# Clear ------
rm(list= ls())
gc()

# Remove the  display limit ------
options(max.print = 100000000)


# Libraries -----
ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}


Packages <- c("tidyverse", "ggplot2", "readstata13", "foreign",  "nnet", "DescTools", 
              "effects", "data.table", "arm", "gdata", "haven", "psy", "stargazer", "mlogit", "eurostat", "survey",
              "rvest", "RStata", "zoo", "padr", "lavaan", "psych", "Hmisc", "labelled", "texreg", "gridExtra")
ipak(Packages)


# Directories -----

# Danilo Mac


load("~/Dropbox/Analisi/dataset_paper_RISP.Rdata")
colnames (dataset_paper)

# Table 1. Multinomial models predicting the intention to vote for the main Italian parties (ref. category Partito Democratico) ------

model <- multinom (Intentionvote_3 ~ Coronavirusknowledge_3 + Vaccine_3 + Trustvaccine_3 +
                     Infosocialmedia_3 + Trstnsp_3 + Trstparl_3 +
                     Lrcat_3 + Performanceeconomy_3 + Economicuncertainty_3 + 
                     Gender_3 + Age_3 + Edulevel_3 + area_3 + Workingstatus_3,  data = dataset_paper, model = TRUE)
screenreg(list(model))


# Table 2. Logistic regression models predicting the intention to vote for populist government and opposition parties (as opposed to mainstream parties) ------
PopOppo<- glm (PopOppo_3 ~ Coronavirusknowledge_3 + Vaccine_3 + Trustvaccine_3 +
                  Infosocialmedia_3 + Trstnsp_3 + Trstparl_3 +
                  Lrcat_3 + Performanceeconomy_3 + Economicuncertainty_3 + 
                  Gender_3 + Age_3 + Edulevel_3 + area_3 + Workingstatus_3,  data = dataset_paper, family = "binomial")


stargazer(PopOppo, type = "text", digits = 2)


PopGov <- glm (PopGov_3 ~ Coronavirusknowledge_3 + Vaccine_3 + Trustvaccine_3 +
                  Infosocialmedia_3 + Trstnsp_3 + Trstparl_3 +
                  Lrcat_3 + Performanceeconomy_3 + Economicuncertainty_3 + 
                  Gender_3 + Age_3 + Edulevel_3 + area_3 + Workingstatus_3,  data = dataset_paper, family = "binomial")

stargazer(PopGov, type = "text", digits = 2)




# Figure 2. Percentage distribution of intention to vote -----
Intention_vote<- dataset_paper %>% 
  drop_na(Intentionvote_3) %>%
  ggplot(aes(x = Intentionvote_3))+
  labs(y="",
       x="",
       title = "Figure 2. Percentage distribution of intention to vote") +
  scale_y_continuous(breaks = seq(0, .25, .05),
                     limits = c(0, .25),
                     labels = scales::percent) +
  geom_text(aes(y = ((..count..)/sum(..count..)), label = scales::percent((..count..)/sum(..count..))), stat = "count", vjust = -0.3) +
  geom_bar(color="white", fill="grey", position = "dodge", aes(y=..count../sum(..count..)))+
  theme_bw() 
Intention_vote


Intention_vote<- Intention_vote + theme(plot.title = element_text (size = 16, face = "bold"),
                                        axis.text.x = element_text(color = "black", size = 9, angle = 0, face = "plain"),
                                        axis.text.y = element_text(color = "black", size = 10, angle = 0, face = "plain"),  
                                        axis.title.x = element_text(color = "black", size = 9, angle = 0, face = "plain"),
                                        axis.title.y = element_text(color = "black", size = 10, angle = 0, face = "plain")) 
Intention_vote



# Figure 4: Coefficient and standard errors plot of the main covariates of interest  ------

Coefficients <- broom::tidy(model,conf.int=TRUE)
Coefficients <- dplyr::filter(Coefficients, term!="(Intercept)")
Coefficients <- dplyr::filter(Coefficients, term!="Trstparl_3")
Coefficients <- dplyr::filter(Coefficients, term!="Infosocialmedia_3")
Coefficients <- dplyr::filter(Coefficients, term!="Trstnsp_3")
Coefficients <- dplyr::filter(Coefficients, term!="Lrcat_3Far left")
Coefficients <- dplyr::filter(Coefficients, term!="Lrcat_3Center-left")
Coefficients <- dplyr::filter(Coefficients, term!="Lrcat_3Center-right")
Coefficients <- dplyr::filter(Coefficients, term!="Lrcat_3Far right")
Coefficients <- dplyr::filter(Coefficients, term!="Lrcat_3I don't identify myself")
Coefficients <- dplyr::filter(Coefficients, term!="Performanceeconomy_3")
Coefficients <- dplyr::filter(Coefficients, term!="Economicuncertainty_3")
Coefficients <- dplyr::filter(Coefficients, term!="Gender_3Female")
Coefficients <- dplyr::filter(Coefficients, term!="Age_3")
Coefficients <- dplyr::filter(Coefficients, term!="Edulevel_3")
Coefficients <- dplyr::filter(Coefficients, term!="area_3North-East")
Coefficients <- dplyr::filter(Coefficients, term!="area_3Center")
Coefficients <- dplyr::filter(Coefficients, term!="area_3South")
Coefficients <- dplyr::filter(Coefficients, term!="area_3Islands")
Coefficients <- dplyr::filter(Coefficients, term!="Workingstatus_3Employed")
Coefficients <- dplyr::filter(Coefficients, term!="area_3North-East")
Coefficients <- dplyr::filter(Coefficients, term!="area_3North-East")

Lega <- dplyr::filter(Coefficients, y.level=="Lega")
M5S <- dplyr::filter(Coefficients, y.level=="Movimento 5 Stelle")
FdI <- dplyr::filter(Coefficients, y.level=="Fratelli d'Italia")

# * Lega ---------------------------------
Plot_Lega<- ggplot(Lega, aes(x=term, y=estimate, ymin=conf.low, ymax=conf.high)) + 
  geom_pointrange() + 
  scale_y_continuous("", breaks=seq(-2, 2, 0.5))  +
  ggtitle("Lega") +
  scale_x_discrete("",
                   limits=c("Coronavirusknowledge_3", "Vaccine_3", "Trustvaccine_3"),
                   labels = c("Index of knowledge about the Coronavirus", "Likelihood to get the vaccination", "Trust in the reliability of the vaccines")) +
  geom_hline(yintercept=0, linetype=2)+
  coord_flip() + theme_classic()

Plot_Lega<- Plot_Lega + theme(plot.title = element_text (size = 16, face = "bold"),
                              axis.text.x = element_text(color = "black", size = 10, angle = 0, hjust = 0, vjust = .5, face = "plain"),
                              axis.text.y = element_text(color = "black", size = 10, angle = 0, hjust = 1, vjust = 0, face = "plain"),  
                              axis.title.x = element_text(color = "black", size = 10, angle = 0, hjust = .5, vjust = .5, face = "plain"),
                              axis.title.y = element_text(color = "black", size = 10, angle = 0, hjust = .5, vjust = .5, face = "plain")) 
Plot_Lega
# * M5S ---------------------------------
Plot_M5S<- ggplot(M5S, aes(x=term, y=estimate, ymin=conf.low, ymax=conf.high)) + 
  geom_pointrange() + 
  scale_y_continuous("", breaks=seq(-2, 2, 0.5))  +
  ggtitle("Movimento 5 Stelle") +
  scale_x_discrete("",
                   limits=c("Coronavirusknowledge_3", "Vaccine_3", "Trustvaccine_3"),
                   labels = c("Index of knowledge about the Coronavirus", "Likelihood to get the vaccination", "Trust in the reliability of the vaccines")) +
  geom_hline(yintercept=0, linetype=2)+
  coord_flip() + theme_classic()

Plot_M5S<- Plot_M5S + theme(plot.title = element_text (size = 16, face = "bold"),
                            axis.text.x = element_text(color = "black", size = 10, angle = 0, hjust = 0, vjust = .5, face = "plain"),
                            axis.text.y = element_text(color = "black", size = 10, angle = 0, hjust = 1, vjust = 0, face = "plain"),  
                            axis.title.x = element_text(color = "black", size = 10, angle = 0, hjust = .5, vjust = .5, face = "plain"),
                            axis.title.y = element_text(color = "black", size = 10, angle = 0, hjust = .5, vjust = .5, face = "plain")) 
Plot_M5S
# * FdI ---------------------------------
Plot_FdI<- ggplot(FdI, aes(x=term, y=estimate, ymin=conf.low, ymax=conf.high)) + 
  geom_pointrange() + 
  scale_y_continuous("", breaks=seq(-2, 2, 0.5))  +
  ggtitle("Fratelli d'Italia") +
  scale_x_discrete("",
                   limits=c("Coronavirusknowledge_3", "Vaccine_3", "Trustvaccine_3"),
                   labels = c("Index of knowledge about the Coronavirus", "Likelihood to get the vaccination", "Trust in the reliability of the vaccines")) +
  geom_hline(yintercept=0, linetype=2)+
  coord_flip() + theme_classic()

Plot_FdI<- Plot_FdI + theme(plot.title = element_text (size = 16, face = "bold"),
                            axis.text.x = element_text(color = "black", size = 10, angle = 0, hjust = 0, vjust = .5, face = "plain"),
                            axis.text.y = element_text(color = "black", size = 10, angle = 0, hjust = 1, vjust = 0, face = "plain"),  
                            axis.title.x = element_text(color = "black", size = 10, angle = 0, hjust = .5, vjust = .5, face = "plain"),
                            axis.title.y = element_text(color = "black", size = 10, angle = 0, hjust = .5, vjust = .5, face = "plain")) 
Plot_FdI

# Aggregate Graphs -----
challenger<-grid.arrange(Plot_Lega, Plot_FdI, Plot_M5S)



# Appendix C: Weighted mean and descriptive statistics  ------
colnames(dataset_paper)
library(pastecs)
stat.desc(dataset_paper)
Hmisc::describe (dataset_paper$Intentionvote_3)
wtd.mean(dataset_paper$Coronavirusknowledge_3, dataset_paper$W3_peso_1, na.rm = TRUE)   
Hmisc::describe (dataset_paper$Coronavirusknowledge_3)
sd (dataset_paper$Coronavirusknowledge_3)
wtd.mean(dataset_paper$Vaccine_3, dataset_paper$W3_peso_1, na.rm = TRUE)   
Hmisc::describe (dataset_paper$Vaccine_3)
sd (dataset_paper$Vaccine_3)
wtd.mean(dataset_paper$Trustvaccine_3, dataset_paper$W3_peso_1, na.rm = TRUE)   
Hmisc::describe (dataset_paper$Trustvaccine_3)
sd (dataset_paper$Trustvaccine_3)

wtd.mean(dataset_paper$Trustvaccine_3, dataset_paper$W3_peso_1, na.rm = TRUE)   
Hmisc::describe (dataset_paper$Infosocialmedia_3)
wtd.mean(dataset_paper$Infosocialmedia_3, dataset_paper$W3_peso_1, na.rm = TRUE)   
Hmisc::describe (dataset_paper$Trstnsp_3)
wtd.mean(dataset_paper$Trstnsp_3, dataset_paper$W3_peso_1, na.rm = TRUE)   
Hmisc::describe (dataset_paper$Trstparl_3)
wtd.mean(dataset_paper$Trstparl_3, dataset_paper$W3_peso_1, na.rm = TRUE)   
Hmisc::describe (dataset_paper$Lrcat_3)
wtd.mean(dataset_paper$Lrcat_3, dataset_paper$W3_peso_1, na.rm = TRUE)   
Hmisc::describe (dataset_paper$Performanceeconomy_3)
wtd.mean(dataset_paper$Performanceeconomy_3, dataset_paper$W3_peso_1, na.rm = TRUE)   
Hmisc::describe (dataset_paper$Economicuncertainty_3)
wtd.mean(dataset_paper$Economicuncertainty_3, dataset_paper$W3_peso_1, na.rm = TRUE)   
Hmisc::describe (dataset_paper$Gender_3)
wtd.mean(dataset_paper$Gender_3, dataset_paper$W3_peso_1, na.rm = TRUE)   
Hmisc::describe (dataset_paper$Age_3)
wtd.mean(dataset_paper$Age_3, dataset_paper$W3_peso_1, na.rm = TRUE)   
Hmisc::describe (dataset_paper$Edulevel_3)
wtd.mean(dataset_paper$Edulevel_3, dataset_paper$W3_peso_1, na.rm = TRUE)   
Hmisc::describe (dataset_paper$area_3)
wtd.mean(dataset_paper$area_3, dataset_paper$W3_peso_1, na.rm = TRUE)   
Hmisc::describe (dataset_paper$Workingstatus_3)
wtd.mean(dataset_paper$Workingstatus_3, dataset_paper$W3_peso_1, na.rm = TRUE)   

Hmisc::describe (dataset_paper$PopOppo_3)
Hmisc::describe (dataset_paper$PopGov_3)



